Análise de Tabelas de contingência
Paulo S. P Silveira (paulo.silveira@fm.usp.br)
José O. Siqueira (jose.siqueira@fm.usp.br)
v20200426.2342
Objetivos
- Identificar e sumariar dados categóricos em tabela de contingência.
- Identificar distribuição amostral \(\chi^2\) (qui-quadrado).
- Formular hipótese estatística na análise de dados categóricos.
- Definir e determinar frequências observada e esperada.
- Identificar, aplicar e interpretar os testes de aderência.
- Identificar e formular análise de variáveis categóricas em tabela de contingência \(l\)x\(c\) (números de linhas versus de colunas).
- Analisar tabelas de contingencia com variáveis ordinais.
- Analisar de tabelas “tridimensionais” \(l\)x\(c\)x\(k\), com a inclusão de variável de estratificação.
- Realizar os procedimentos estatísticos em R.
Preparação
Para executar os exemplos aqui apresentados, sugerimos que crie um projeto e baixe os seguintes arquivos para a pasta do mesmo:
- tabelacontingencia2x2_capacetetrauma.R
- Qui2.R
- Qui2_LxC.R
- friendlycolor.R
- TesteQuiQuadradoAderencia_Chocolate.R
- chocolate.xlsx
- eiras.quiaderencia.R
- TesteQuiQuadradoAderencia_Flores.R
- flores.xlsx
- TabelaContingencia2x2_TabagismoEtilismo.R
- tabagismo_e_etilismo_2x2.xlsx
- tabagismo_e_etilismo_2x2_v2.xlsx
- TCrXc_StatusOcupacional.R
- Teste qui-quadrado de Cochran-Armitage.R
- menarca.xlsx
- Teste qui-quadrado de Cochran-Armitage generalizado.R
- Teste Qui-quadrado de Mantel-Haenszel.R
- fumo_e_faixaetaria.xlsx
Roteiro
Há diversas funções e procedimentos envolvidos neste texto. Como roteiro orientador, verifique:
Delineamento entre participantes
teste de independência / dissociação
- Tabela unidimensional
- Uma variável nominal:
- Teste de aderência: chisq.test
- Uma variável nominal:
- Tabela bidimensional
- Duas variáveis nominais dicotômicas
- Teste qui-quadrado de Pearson: chisq.test
- Teste exato de Fisher de OR bruto: fisher.test
- Teste qui-quadrado mínimo
- Uma variável nominal dicotômica e uma ordinal
- Teste qui-quadrado de Cochran-Armitage (trend): - - DescTools::CochranArmitageTest
- Uma variável nominal politômica e uma ordinal
- Teste qui-quadrado de Cochran-Armitage generalizado: coin::chisq_test
- Duas variáveis ordinais
- Teste qui-quadrado de Mantel-Haenszel: DescTools::MHChisqTest
- Duas variáveis nominais dicotômicas
- Tabela tridimensional
- Duas variáveis nominais dicotômicas e uma variável de estratificação nominal politômica
- Teste de Mantel-Haenszel de OR comum: stats::mantelhaen.test
- Duas variáveis nominais politômicas e uma variável de estratificação nominal politômica
- Teste de Mantel-Haenszel de OR comum generalizado: coin::cmh_test
- Duas variáveis ordinais e uma variável de estratificação nominal politômica
- Teste qui-quadrado ordinal estratificado (linear-by-linear association model): coin::lbl_test
- Duas variáveis nominais dicotômicas e uma variável de estratificação nominal politômica
Delineamento intraparticipantes
teste de concordância entre dois métodos de avaliação
- Tabela bidimensional
- Duas variáveis nominais dicotômicas
- Teste de McNemar: mcnemar.test, coin::mh_test
- Teste kappa de Cohen: epiR::epi.kappa
- Duas variáveis nominais dicotômicas
Contingência
Em sentido geral, contingência pode significar qualquer relação de dependência entre eventos ambientais ou entre eventos comportamentais e ambientais (Catania, 1993; Skinner, 1953, 1969; Todorov, 1985). Embora possa ser encontrado nos dicionários com diferentes significados, esse termo é empregado, na análise do comportamento, como termo técnico para enfatizar como a probabilidade de um evento pode ser afetada ou causada por outros eventos (Catania, 1993, p. 368).http://www.itcrcampinas.com.br/txt/texto_deisy.pdf
\[~\]
Relações entre duas variáveis
|
1 Exposição e Desfecho |
tabela 2x2
independência
concordância
discordância
\[~\]
Associação
|
2 Não avalia causa e efeito: apenas associação. |
causa e efeito positivo?
causa e efeito negativo?
causa e efeito positivo?
Exemplos
Gatos e alergia
- estresse
- mecanismo imune
- desinfetantes
- alergia e estresse
- alergia e gatos
Vinho e doença coronariana
- ambiente inseguro
- moradia precária
- má nutrição
- boa moradia
- segurança
- boa alimentação
- outros cuidados
- é o vinho?
Fumo e Câncer de Pulmão
- The Cigarette Century
- How could one attribute lung cancer to cigarette smoking?
“The public impression of scientific and medical uncertainty that resulted was a crucial element in maintaining the market of current smokers as well as recruiting new ones. Industry literature, for example, frequently pointed to the fact that nonsmokers could also develop lung cancer. Therefore, they argued, how could one attribute lung cancer to cigarette smoking?”
- Medical knowledge is provisional
medical knowledge was always provisional and contingent. Just as drugs deemed “effective” do not work in every case, so too a cause of disease does not always result in disease. As Richard Doll would later explain, the epidemiologists had identified a cause of lung cancer (and other diseases), not the cause.
- Fisher was against the association between smoking and lung cancer
Another vocal critic of the lung cancer findings was Sir Ronald Fisher, the leading biometrician and geneticist in Great Britain during the first half of the twentieth century and a man deeply committed to bringing statistical analysis to genetics and agricultural experimentation. His 1925 book, Statistical Methods for Research Workers, quickly became a classic, leading to academic appointments at University College London and Cambridge University. Fisher’s critiques were similar to Berkson’s. The ethical impossibility of conducting a randomized experiment led him to question the results of the epidemiological studies. As a believer in genetic* notions of cancer causality, Fisher speculated that there was some constitutional factor that led individuals both to become smokers and to get lung cancer, even though smoking and lung cancer might not be causally related. Doll and Hill repeatedly rebutted this theory, returning to the critical question of how to account for the rise in lung cancers during the twentieth century if the disease was simply “constitutional.”
* 1923: One of the [A. Bradford] Hill’s innovations was the first randomized, double-blind clinical trial, designed to reduce investigator bias in the evaluation of clinical outcomes. […] This method, which drew on Fisher’s agricultural experimentation in genetics, became a critical new tool for evaluating medical treatments.
By the 1930s, the relationship of smoking to cancer was a topic of uresolved debate. It was the life insurance industry, which like the tobacco corporations grew by leaps and bounds in the first half of the twentieth century, that took the lead in understanding the effects of smoking on health.
By the early 1950s, however, it was abundantly clear that the evidence implicating cigarette smoking as a risk to health was now of a different order. First, the link between smoking and disease was categorical, outside the realm of individual clinical judgment.
They [tobacco industry] responded with a new and unprecedent public relations strategy. Its goals was to produce and sustain scientific skepticism and controversy in order to disrupt the emerging consensus on the harms of cigarette smoking. This strategy required intrusions into scientific process and procedure. […] So long as there appeared to be doubt, so long as the industry could assert “not proven”. […] The future of the cigarette would now depend on the successful production of a scientific controversy.
Óculos e xixi
Quando alguém lhe traz uma opinião coincidente com sua crença, você a toma como um fato. Quando alguém lhe traz um fato discordante de sua crença, você o toma como uma opinião.
A hipótese nula
|
3 Tomada de decisão: Hipótese nula (\(H_0\)) e hipótese alternativa (\(H_1\)) |
|
Adotaremos, neste texto, a seguinte convenção:
|
Teste Qui-quadrado
Os testes baseados na estatística qui-quadrado podem ser divididos em dois tipos:
Teste de aderência uma variável nominal com duas ou mais categorias: teste de frequências hipotetizadas
Teste de independência entre duas variáveis nominais com duas ou mais categorias
- a distribuição \(\chi^2\)
Antes de discutirmos os testes, vamos explorar um pouco a distribuição que utilizam. Diferente da distribuição normal que têm domínio de \([- \infty, + \infty]\), estas são distribuições assimétricas, iniciadas em zero, que tendem assintoticamente a zero com o valor crescente de \(\chi^2\), \([0, + \infty]\).
O comando para calcular a probabilidade acumulada entre zero e cada valor de qui-quadrado é:
dchisq(x, df, ncp = 0, log = FALSE)
onde x é o valor de \(\chi^2\), df são os graus de liberdade (degrees of freedom) e ncp é o parâmetro de não centralidade (non-centrality parameter = 0, por default, correspondente à distribuição centrada de \(H_0\); voltaremos ao ncp adiante). O parâmetro log, desligado por default é para a possibilidade de devolver o logarítmo do valor (não o exploraremos neste texto).
Pode-se gerar curvas desta distribuição com o seguinte código, mostrando seu aspecto entre 2 e 6 graus de liberdade:
chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 2:6
for (g in 1:length(graus.liberdade))
{
p <- dchisq(chi2.valor, df=graus.liberdade[g])
plot(chi2.valor, p, type="l",
main=paste("Distribuição de Qui^2, df=",
graus.liberdade[g],sep=""),
xlab="Qui-quadrado", ylab="Densidade")
}
ou, caso prefira ver todas as curvas no mesmo gráfico, sugere-se:
source("friendlycolor.R")
cor <- 3
padrao <- 1
leg.txt <- c()
leg.cor <- c()
chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 2:6
for (g in 1:length(graus.liberdade))
{
p <- dchisq(chi2.valor, df=graus.liberdade[g])
if (g==1)
{
plot(chi2.valor, p, type="l", lty=padrao,
col=friendlycolor(cor), lwd=2,
main="Distribuições de Qui^2",
xlab="Qui-quadrado", ylab="Densidade")
} else
{
lines(chi2.valor, p, lty=padrao,
col=friendlycolor(cor), lwd=2)
}
# guarda para a legenda
leg.txt <- c(leg.txt, paste("df = ",graus.liberdade[g],sep=""))
leg.cor <- c(leg.cor, friendlycolor(cor))
cor <- cor+6
padrao <- padrao+1
}
legend("topright",
leg.txt,
col=leg.cor,
lwd=2,
lty=1:5,
box.lwd=0, bg="transparent") 
Teste de aderência: teste qui-quadrado de uma variável
Permite descobrir se um conjunto de frequências observadas difere de um outro conjunto de frequências hipotetizadas.
suposição
Independência das observações: cada observação da unidade observacional deve ser alocada em apenas um categoria da variável nominal.exemplo 1: chocolate
Cento e dez pessoas foram solicitadas a manifestar suas preferências com respeito a chocolates. Queremos verificar se alguma das marcas é preferida em detrimento de outras.
Se não são (\(H_0\)), devemos esperar aproximadamente o mesmo número de pessoas em cada categoria.
Com certeza, não haverá exatamente o mesmo número de participantes em cada categoria, mas os números devem ter uma distribuição com razoável proximidade.
O número de pessoas que escolheu apenas uma entre quatro marcas diferentes está em chocolate.xlsx.
O teste está implementado em TesteQuiQuadradoAderencia_Chocolate.R:
# TesteQuiQuadradoAderencia_Chocolate.R
library(readxl)
cat("Teste de aderencia a distribuicao multinomial\n")
df_chocolate <- read_excel("chocolate.xlsx")
prmatrix(df_chocolate, rowlab = rep("",nrow(df_chocolate)))
out <- chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0,
simulate.p.value = TRUE, B = 1e6)
df <- chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0)$parameter
X2 <- out$statistic
p <- out$p.value
n <- sum(chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0)$observed)
V <- sqrt((X2/n)/df)
cat("X^2 =", X2 , "gl =", df, "p exato =", p, "\n")
cat("\nHeuristica de significancia: X^2/gl > 2:\n", X2/df, "\n")
cat("\nResiduos estandardizados:\n", out$stdres, "\n")
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("V de Cramer = ", V,"\n",sep="")
cat("(Grau ", gV, "de dessemelhanca entre as distribuicoes observada e hipotetizada\n",sep="")Observe o código:
- lê a planilha no formato xlsx com os dados e a exibe.
- note que denominamos uma coluna como probH0, contendo \(0.25\) em todas as linhas; estas são as frequências esperadas para a hipótese nula, \(H_0\).
- O teste em si recebe a coluna com o número de pessoas que preferiram cada chocolate e estas frequências esperadas; é executado com bootstrapping, com \(10^5\) iterações (parâmetro \(B\)).
- O restante do código serve para reorganizar e ecoar os resultados, incluindo o cálculo de tamanho de efeito dado pelo V de Cramér.
Este código gera a seguinte saída:
Teste de aderencia a distribuicao multinomial
Chocolate Pessoas probH0
"A" "20" "0.25"
"B" "60" "0.25"
"C" "10" "0.25"
"D" "20" "0.25"
X^2 = 53.63636 gl = 3 p exato = 9.99999e-07
Heuristica de significancia: X^2/gl > 2:
17.87879
Residuos estandardizados:
-1.651446 7.156264 -3.853373 -1.651446
V de Cramer = 0.4031556
(Grau intermediariode dessemelhanca entre as distribuicoes observada e hipotetizadaCom \(p=9.999 \cdot 10^-5\), rejeita-se \(H_0\) e, portanto, temos evidência de que a escolha das pessoas pelas marcas de chocolate não é homogênea. Observando a tabela, vemos que o desequilíbrio ocorreu pela maior preferência pelo chocolate B.
| Os graus de liberdade em uma tabela com \(L\) linhas e \(C\) colunas são dados por \[df = (L-1) \cdot (C-1)\] Portanto, uma tabela \(2\text{ x }2\), tem somente um grau de liberdade, e uma tabela \(4\text{ x }3\) tem seis graus de liberdade. |
exemplo 2: genética
As probabilidades relacionadas com \(H_0\) não precisam ser sempre as mesmas. Por exemplo, em genética, é comum precisarmos avaliar como é a herança de um genótipo a partir dos fenótipos observados.
Duas linhagens puras de plantas foram cruzadas: uma com pétalas amarelas e outra com pétalas vermelhas. A primeira geração, \(F_1\), tem pétalas cor de laranja. Então \(F_1\) é cruzada entre si, gerando \(F_2\) com 320 plantas, das quais 182 têm pétalas cor de laranja, 61 amarelas e 77 vermelhas.Há duas possibilidades para explicar este resultado, e queremos saber qual é o caso para o fenótipo das pétalas das flores desta espécie:
- dominância incompleta ou codominância, com os seguintes genótipos e frequências esperadas:
- Y/R: laranja (\(1 \over 2\))
- Y/Y: amarela (\(1 \over 4\))
- R/R: vermelha (\(1 \over 4\))
Neste caso, os números esperados para \(F_2\) são, respectivamente, 160, 80 e 80, seguindo as proporções 2:1:1.
- epistasis recessiva, com os seguintes genótipos e frequências esperadas:
- Y/-; R/-: laranja (\(9 \over 16\))
- y/y; R/-: amarela (\(3 \over 16\))
- Y/-; r/r: vermelha (\(3 \over 16\))
- y/y; r/r: vermelha (\(1 \over 16\))
Neste caso, os números esperados para \(F_2\) são, respectivamente, 180 , 60 e 80 (60+20), seguindo as proporções 9:3:(3:1).
A qual das duas hipóteses a herança da cor das flores desta espécie adere?
Implementamos uma solução em eiras.quiaderencia.R e TesteQuiQuadradoAderencia_Flores.R.
Observe o código com nova função declarada:
# eiras.quiaderencia.R
# Teste Qui-quadrado de aderencia
eiras.quiaderencia <- function (df_dados, B=1e6)
{
valores <- as.vector(unlist(df_dados[,1]))
probabilidades <- as.vector(unlist(df_dados[,2]))
# teste robusto (bootstrapping)
robusto <- chisq.test(valores, p=probabilidades,
simulate.p.value = TRUE, B = B)
# usando o teste convencional para obter os graus de liberdade
tradicional <- chisq.test(valores, p=probabilidades)
df <- tradicional$parameter
X2 <- robusto$statistic
p <- robusto$p.value
n <- sum(robusto$observed)
V <- sqrt((X2/n)/df)
# resultados
cat("----------------------------------------------------\n")
cat("Teste de aderencia a distribuicao multinomial\n")
prmatrix(df_dados, rowlab = rep("",nrow(df_dados)))
cat("X^2 =", X2 , "gl =", df, "p exato =", p, "\n")
cat("\nHeuristica de significancia: X^2/gl > 2:\n", X2/df, "\n")
cat("\nResiduos estandardizados:\n", robusto$stdres, "\n")
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("V de Cramer = ", V,"\n",sep="")
cat("(Grau ", gV, " de dessemelhanca entre as distribuicoes observada e hipotetizada\n",sep="")
cat("----------------------------------------------------\n")
}Criar uma função e armazená-la em um arquivo separado facilita seu uso (o nome do arquivo tem o mesmo nome da função só por conveniência; poderia ter qualquer nome). Esta função, eiras.quiaderencia(), recebe um data frame com duas colunas e o número de iterações para o bootstrapping (com \(10^6\) por default).
Adaptar o código para outros casos torna-se muito mais simples. No código anterior, seria preciso alterar todas as ocorrências de df_chocolate e os nomes das colunas Pessoas e probH0. Agora, só é necessário modificar TesteQuiQuadradoAderencia_Flores.R:
# TesteQuiQuadradoAderencia_Flores.R
library(readxl)
source ("eiras.quiaderencia.R")
# pega o arquivo e mostra seu conteúdo
df_flores <- read_excel("flores.xlsx")
cat("\nDados:\n")
prmatrix(df_flores, rowlab = rep("",nrow(df_flores)))
# chama os testes estatísticos
cat("\nHipotese de codominancia:\n")
eiras.quiaderencia(df_flores[,c(2,3)], B=1e6)
cat("\nHipotese de epistasis recessiva:\n")
eiras.quiaderencia(df_flores[,c(2,4)], B=1e6)Para funcionar, carrega-se a função na memória (com source()), lemos o arquivo de interesse flores.xlsx e chamamos eiras.quiaderencia() duas vezes, indicando as colunas a serem testadas:
- a coluna 2 da planilha (que tem o número de plantas com cada característica) contra a coluna 3, que testa a hipótese de codominância,
- a coluna 2 da planilha contra a coluna 4, que testa a hipótese de epistasis.
Este código produz a seguinte saída:
Dados:
CorPetalas Numero Codominancia Epistasis
"laranja" "182" "0.50" "0.5625"
"vermelha" " 77" "0.25" "0.2500"
"amarela" " 61" "0.25" "0.1875"
Hipotese de codominancia:
----------------------------------------------------
Teste de aderencia a distribuicao multinomial
Numero Codominancia
182 0.50
77 0.25
61 0.25
X^2 = 7.65 gl = 2 p exato = 0.02224998
Heuristica de significancia: X^2/gl > 2:
3.825
Residuos estandardizados:
2.459675 -0.3872983 -2.452889
V de Cramer = 0.1093303
(Grau pequeno de dessemelhanca entre as distribuicoes observada e hipotetizada
----------------------------------------------------
Hipotese de epistasis recessiva:
----------------------------------------------------
Teste de aderencia a distribuicao multinomial
Numero Epistasis
182 0.5625
77 0.2500
61 0.1875
X^2 = 0.1513889 gl = 2 p exato = 0.9324391
Heuristica de significancia: X^2/gl > 2:
0.07569444
Residuos estandardizados:
0.2253745 -0.3872983 0.143223
V de Cramer = 0.01538002
(Grau minimo de dessemelhanca entre as distribuicoes observada e hipotetizada
----------------------------------------------------Neste exemplo buscamos não rejeitar \(H_0\). Para o habitual \(\alpha=0.05\) estes resultados sugerem que a herança das cores das pétalas desta espécie é explicada por epistasis recessiva.
Teste Qui-quadrado de Pearson para Independência
O teste Qui-quadrado de Pearson permite que se descubra se existe um relacionamento ou efeito de interação entre duas variáveis qualitativas nominais com duas ou mais categorias.
Cada observação da unidade observacional deve ser alocada em apenas uma categoria de cada uma das duas variáveis nominais.
Exemplo: capacete e trauma
Considere um estudo que investiga a efetividade dos capacetes de segurança de bicicleta na prevenção de traumas cranianos. Os dados consistem de uma amostra aleatória de 793 indivíduos envolvidos em acidentes ciclísticos durante um certo período. Deseja-se testar, com \(\alpha=0.05\), se o uso do capacete tem funcionado como um fator de proteção.
Formulamos as hipóteses nula e alternativa:
\(H_0: \text{não existe associação entre uso do capacete e trauma craniano}\)
\(H_1: \text{existe associação entre uso do capacete e trauma craniano}\)
\(\alpha = 5\%\)
Os dados são: 
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
print (tabela) Trauma + Trauma -
Capacete + 17 138
Capacete - 130 508Sob a hipótese nula (\(H_0\)), i.e., se não existe associação entre o uso de capacete e o trauma craniano nos acidentes, os valores esperados são:
chi2 <- chisq.test(tabela)
print(round(chi2$expected,3)) Trauma + Trauma -
Capacete + 28.733 126.267
Capacete - 118.267 519.733É possível transformar em proporções, obtendo:
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
tabela_prop <- tabela/sum(tabela)
print (round(tabela_prop,3)) Trauma + Trauma -
Capacete + 0.021 0.174
Capacete - 0.164 0.641e, sob \(H_0\), os valores esperados são:
soma_col1 <- sum(tabela_prop[,1])
soma_col2 <- sum(tabela_prop[,2])
soma_lin1 <- sum(tabela_prop[1,])
soma_lin2 <- sum(tabela_prop[2,])
tabela_propesp <- tabela_prop
tabela_propesp[1,1] <- soma_lin1*soma_col1
tabela_propesp[1,2] <- soma_lin1*soma_col2
tabela_propesp[2,1] <- soma_lin2*soma_col1
tabela_propesp[2,2] <- soma_lin2*soma_col2
print(round(tabela_propesp,3)) Trauma + Trauma -
Capacete + 0.036 0.159
Capacete - 0.149 0.655Note que o teste em R é feito diretamente com os valores observados (pois o tamanho da amostra importa) e seus resultados foram armazenados, acima, na variável chi2. O tipo desta variável é dado por:
str(chi2)List of 9
$ statistic: Named num 6.7
..- attr(*, "names")= chr "X-squared"
$ parameter: Named int 1
..- attr(*, "names")= chr "df"
$ p.value : num 0.00964
$ method : chr "Pearson's Chi-squared test with Yates' continuity correction"
$ data.name: chr "tabela"
$ observed : 'table' num [1:2, 1:2] 17 130 138 508
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ expected : num [1:2, 1:2] 28.7 118.3 126.3 519.7
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ residuals: 'table' num [1:2, 1:2] -2.189 1.079 1.044 -0.515
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ stdres : 'table' num [1:2, 1:2] -2.7 2.7 2.7 -2.7
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
- attr(*, "class")= chr "htest"sapply(chi2, typeof) statistic parameter p.value method data.name observed
"double" "integer" "double" "character" "character" "double"
expected residuals stdres
"double" "double" "double" É uma estrutura complexa, da qual muitos valores podem ser extraídos quando quisermos utilizá-los dentro de um RScript. A parte em chi2$expected, que já utilizamos, contém a tabela com os valores esperados sob \(H_0\).
O resultado principal do teste é exibido pela variável como um todo:
print (chi2)
Pearson's Chi-squared test with Yates' continuity correction
data: tabela
X-squared = 6.7001, df = 1, p-value = 0.009641A estatística de teste tem o valor \(X^2=6.7001\), com 1 grau de liberdade (\(df=1\)), correspondendo ao valor-p de \(0.009641\) (que será discutido adiante). Com estes valores podemos decidir pela rejeição ou não-rejeição de \(H_0\).
No exemplo de uma tabela 2x2 há somente 1 grau de liberdade. O gráfico da distribuição de qui-quadrado (para melhor visibilidade, computamos \(0 \le \chi^2 \le 8\), pois o valor de \(\chi^2\) vai de zero a \(\infty\), e truncamos o eixo \(y\) em 1.2) é:
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
p <- dchisq(chi2.valor, df=graus.liberdade)
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
xlab="Qui-quadrado", ylab="Densidade"
)
|
Conceitualmente, o teste Qui-quadrado é um teste de proporções, no qual a estatística de teste (i.e., o valor de \(X^2\)) é obtida pela diferença quadrática entre as frequências absolutas observadas e esperadas, dada por \[X^2 = \sum{ {(o - e)^2} \over{e} }\]. ou \[X^2 = {{(ad-bd)^2(a+b+c+d)} \over {(a+b)(c+d)(b+d)(a+c)}}\] Em R, experimente calcular a estatística na “força bruta”: O valor de \(X^2\) aumenta com a discrepância crescente entre os valores observados e esperados. Note, porém, que o valor calculado pelo código acima é diferente do que aparece em: Isto acontece por causa da correção para continuidade de Yates, que será comentada adiante. Na “força bruta”, a correção pode ser implementada com: |
Nosso interesse é saber a partir de qual o valor de \(X^2\) a discrepância entre os valores observados e esperados é menos provável que o valor de \(\alpha\) escolhido. Isto é encontrar o valor que separa a área de \(5\%\) sob a curva de distribuição em sua cauda direita. Este valor é chamado de qui-quadrado crítico, \(\chi^2_c\).
| Note que o raciocínio é bicaudal, buscando associação positiva ou negativa, mas a operacionalização do teste só utiliza a cauda superior. Isto acontece porque \(chi^2=0\) indica igualdade numérica (e, portanto, também estatística) entre os valores esperados e observados. Consequentemente, procuramos \(\alpha=0.05\) somente na cauda direita. |
Em R, podemos encontrar \(\chi^2_c\) pelo complemento (\(1-\alpha=0.95\)) com
chi2.critico <- qchisq(p=0.95, df=1)
cat("qui-quadrado crítico (df=1, alfa=5%): ", chi2.critico, "\n", sep="")qui-quadrado crítico (df=1, alfa=5%): 3.841459pois este é o valor de \(\chi^2\) que deixa \(95\%\) da área sob a curva à sua esquerda.
Graficamente, podemos hachurar em azul a área correspondente a \(\alpha\) e traçar uma linha pontilhada vertical em \(\chi^2_c\) com:
source("friendlycolor.R")
# vetor com valores de 0 a 8
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
# vetor com probabilidades correspondentes
p <- dchisq(chi2.valor, df=graus.liberdade)
# exibe o gráfico
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
xlab="Qui-quadrado", ylab="Densidade")
# vetor com os valores acima de qui critico
chi2.critico <- qchisq(p=0.95, df=1)
chi2maioralfa <- chi2.valor[chi2.valor>chi2.critico]
# probabilidades correspondentes
pmaioralfa <- dchisq(chi2maioralfa, df=graus.liberdade)
# hachura sob a curva
chi2maioralfa <- c(min(chi2maioralfa), chi2maioralfa, max(chi2maioralfa))
pmaioralfa <- c(0, pmaioralfa, 0)
polygon(chi2maioralfa, pmaioralfa, col=friendlycolor(10), border = NA)
# reforca a linha
lines(chi2maioralfa[2:(length(chi2maioralfa)-1)],
pmaioralfa[2:(length(pmaioralfa)-1)],
col=friendlycolor(10), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2.critico, lty=2)
Esta área hachurada é a região de rejeição de \(H_0\).
- decisão pela estatística de teste
A estatística \(\chi^2\) é maior quanto maior for a discrepância entre os valores observados e esperados. Para um grau de liberdade e \(\alpha=0.05\), computamos \(\chi^2_c=\) 3.8414588. No exemplo de capacete e trauma, acima, rejeitamos \(H_0\) porque calculamos \(X^2=\) 6.7001206, valor maior que o valor crítico.
- decisão pelo valor-p
À direita do valor calculado \(X^2=\) 6.7001206 define-se uma área, que corresponde ao valor-p. Este valor também foi calculado quando aplicamos o teste e guardamos seu resultado em chi2, acima; o valor-p encontra-se acessível em chi2$p.value=0.0096406. Como \(p < \alpha\), concluimos pela rejeição de \(H_0\).
Podemos hachurar em laranja a área correspondente ao valor-\(p\) e traçar uma linha pontilhada vertical na altura do \(X^2\) calculado, adicione mais linhas ao código anterior, mas precisamos da variável chi2$statistic que retorna de chisq.tes(). O gráfico, agora, tem muitos elementos e, portanto, também devemos adicionar uma legenda.
O código completo é:
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
chi2 <- chisq.test(tabela)source("friendlycolor.R")
# vetor com valores de 0 a 8
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
# vetor com probabilidades correspondentes
p <- dchisq(chi2.valor, df=graus.liberdade)
# exibe o gráfico
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
xlab="Qui-quadrado", ylab="Densidade")
# vetor com os valores acima de qui critico
chi2.critico <- qchisq(p=0.95, df=1)
chi2maioralfa <- chi2.valor[chi2.valor>chi2.critico]
# probabilidades correspondentes
pmaioralfa <- dchisq(chi2maioralfa, df=graus.liberdade)
# hachura sob a curva
chi2maioralfa <- c(min(chi2maioralfa), chi2maioralfa, max(chi2maioralfa))
pmaioralfa <- c(0, pmaioralfa, 0)
polygon(chi2maioralfa, pmaioralfa, col=friendlycolor(10), border = NA)
# reforca a linha
lines(chi2maioralfa[2:(length(chi2maioralfa)-1)],
pmaioralfa[2:(length(pmaioralfa)-1)],
col=friendlycolor(10), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2.critico, lty=2)
# ########## Area a direita de p ############
# vetor com os valores acima de qui calculado
chi2maiorp <- chi2.valor[chi2.valor>chi2$statistic]
# probabilidades correspondentes
pmaiorp <- dchisq(chi2maiorp, df=graus.liberdade)
# hachura sob a curva
chi2maiorp <- c(min(chi2maiorp), chi2maiorp, max(chi2maiorp))
pmaiorp <- c(0, pmaiorp, 0)
polygon(chi2maiorp, pmaiorp, col=friendlycolor(20), border = NA)
# reforca a linha
lines(chi2maiorp[2:(length(chi2maiorp)-1)],
pmaiorp[2:(length(pmaiorp)-1)],
col=friendlycolor(20), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2$statistic, lty=3)
# legenda
legend("topright",
c("Qui^2 crítico","Qui^2 obs.", "alfa", "p"),
col=c("black", "black", friendlycolor(10), friendlycolor(20)),
lty=c(2, 3, 1, 1),
lwd=c(1, 1, 5, 5),
box.lwd=0, bg="white",
cex=0.8) 
|
As duas formas de decisão são igualmente válidas e equivalentes. Exibir o valor-p é a forma mais moderna e recomendada atualmente. |
Teste Qui-quadrado robusto
Repare nas saídas da função chisq.test(), cujo cabeçalho indica:
Pearson’s Chi-squared test with Yates’ continuity correction
Esta é uma versão tradicional do teste, para a qual a correção de Yates é aplicável por default às tabelas 2x2 (mais detalhes serão vistos adiante). Podemos evitá-la com:
chi2ny <- chisq.test(tabela, correct=FALSE)
print(chi2ny)
Pearson's Chi-squared test
data: tabela
X-squared = 7.3099, df = 1, p-value = 0.006858Além disto, há várias condições para a aplicação do teste Qui-quadrado (comentadas adiante) que podem, muitas vezes, restringir seu uso.
Uma alternativa é computar o teste em sua versão robusta. A função chisq.test() consegue simular a distribuição por bootstrapping, alterando-se o comando para:
chi2.robusto <- chisq.test(tabela, correct=FALSE,
simulate.p.value = TRUE, B=1e6)
print(chi2.robusto)
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tabela
X-squared = 7.3099, df = NA, p-value = 0.007638Esta versão é simulada com \(1e6 = 1 \cdot 10^6 = 1000000\) iterações.
| Como a amostra é grande, os valores de \(X^2\) e \(p\) não mudam muito entre as diferentes versões. Há um pequeno bug na implementação da versão robusta, que não calcula os graus de liberdade (sabemos que \(df=1\)). |
A decisão, neste caso, não mudou, tanto pelo critério do \(\chi^2_c\) quanto pelo \(p\): rejeita-se \(H_0\).
Riscos relativos
|
Odds e probabilidade Odds (traduzido habitualmente por chance) e probabilidade são formas equivalentes para expressar possibilidades e relacionadas por: \(\textit{Odds} = {{\text{Probabilidade}} \over{1-\text{Probabilidade}}}\) e \(\text{Probabilidade} = {{\textit{Odds}} \over{1+\textit{Odds}}}\) Por exemplo, uma probabilidade de \(80\%\) corresponde a \(\textit{Odds} = {{0.8} \over {1-0.8}} = {0.8 \over 0.2} = 4\) indicando que dizer “\(80\%\) de probabilidade” é equivalente a dizer “4 vezes mais chance de ocorrer do que não ocorrer” um determinado evento. Resersamente, Odds de 2 é: \(\text{Probabilidade} = {2 \over {1+2}} = {2 \over 3} \approx 66.67\%\) e \(\textit{Odds}=1\) resulta em: \(\text{Probabilidade} = {1 \over {1+1}} = {1 \over 2} = 50\%\) Associa-se o máximo de incerteza com probabilidade de \(50\%\); por este motivo, \(Odds=1\) será um valor necessário para decisões estatísticas, adiante. |
Há dois tipos principais de risco relativo: razão de riscos (RR, risk ratio) e razão de chances (OR, odds ratio).
Vamos considerar a seguinte tabela:
Tem efeito Não tem efeito
Tem exposição a b a+b
Não tem exposição c d c+d
a+c b+d totalEsta tabela de contingência genericamente relaciona exposição (e.g., hábito de fumar) com o efeito (e.g., câncer de pulmão). Traz as contagens dos:
- expostos que apresentaram o efeito (\(a\)),
- expostos que não apresentaram o efeito (\(b\)),
- não expostos que apresentaram o efeito (\(c\)) e
- não expostos que não apresentaram o efeito (\(d\)),
razão de riscos (RR, risk ratio)
Riscos são probabilidades. A razão de riscos é dada por:
\(RR = \LARGE{{{{a} \over {a+b}} \over {{{c} \over {c+d}}}}}\)
Repare que \({a} \over {a+b}\) é a probabilidade de um indivíduo exposto ter o efeito (e.g., fumante ter câncer); \({c} \over {c+d}\) é a probabilidade de um indivíduo não exposto ter o efeito (e.g., não-fumante ter câncer).
RR, portanto, é um odds: quantas vezes mais é esperado o efeito (ocorrência de câncer) de quem é exposto (fumante) em comparação com quem não é exposto (não fumante);
razão de chances (OR, odds ratio)
Odds ratio como diz o nome, é uma razão entre dois odds, dada por:
\(OR = \LARGE{{ { { {a} \over {a+c} } \over { {c} \over {a+c} } } \over { { {b} \over {b+d} } \over { {d} \over {b+d} } } }}\)
Considere, agora, que:
- \(\LARGE{ { {a} \over {a+c} } \over { {c} \over {a+c} } }\), entre os doentes, quantas vezes é a maior a chance (odds) de ter efeito (câncer) para quem é exposto (fuma).
- \(\LARGE{ { {b} \over {b+d} } \over { {d} \over {b+d} } }\), entre os não doentes, quantas vezes é a maior a chance (odds) de ter efeito (câncer) para quem é exposto (fuma).
Interessante, também, é que (com alguma álgebra) \(OR\) reduz-se a:
\(OR = \LARGE{ {{a \cdot d} \over {{b \cdot c}} } }\)
fácil de lembrar e assim conhecido como “produto cruzado” (com o defeito de ocultar o motivo pelo qual é um odds ratio). Por outro lado, como produto cruzado, percebe-se que é uma medida de quanto pesa a:
- diagonal principal, do que é concordante: \(\text{efeito entre expostos}~~\text{E}~~\text{não efeito entre não expostos}\)
em relação à
- diagonal secundária, do que é discordante: \(\text{efeito entre não expostos}~~\text{E}~~\text{não efeito entre expostos}\).
| Os valores de \(RR\) e \(OR\) são próximos quando menos que \(5\%\) dos não expostos têm efeito. |
Em R, o cálculo de \(OR\) aparece no Teste Exato de Fisher. No exemplo de capacete e trauma, temos:
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
print(tabela)
ft <- fisher.test(tabela) # Teste de OR robusto
print(ft) Trauma + Trauma -
Capacete + 17 138
Capacete - 130 508
Fisher's Exact Test for Count Data
data: tabela
p-value = 0.005688
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.2630993 0.8357342
sample estimates:
odds ratio
0.4817645 - decisão pelo intervalo de confiança
A estimativa pontual de odds ratio está em ft$estimate=0.4817645, mas para a decisão é obrigatório verificar seu intervalo de confiança 95% (\(IC95\)), com limites inferior e superior respectivamente guardados em ft$conf.int[1]=0.2630993 e ft$conf.int[2]=0.8357342.
Como a ausência de efeito é dada por \(OR=1\) (correspondendo a \(H_0\)) e este valor unitário está fora do intervalo, rejeita-se \(H_0\). Além disto, o \(IC95\) é menor que o esperado por \(H_0\) e, portanto, o uso de capacete está associado com redução dos traumas cranianos.
|
Caso a tabela fosse construída com as colunas em posições trocadas, obter-se-ia o seguinte: O odds ratio é agora maior que o valor unitário, o qual está fora do intervalo, rejeitando-se a hipótese nula e a associação é positiva: a conclusão é a mesma, o uso de capacete está associado com o aumento de “Trauma -”, portanto é protetor. Aliás, observe que |
\(RR\) e \(OR\) em epitools
Há outros pacotes em R que calculam risk ratio e odds ratio, como por exemplo:
library(epitools)
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
epitools::oddsratio(tabela, method="fisher")que produz os mesmos valores para a estimativa pontual e intervalo de confiança que o teste exato de Fisher:
$data
Trauma + Trauma - Total
Capacete + 17 138 155
Capacete - 130 508 638
Total 147 646 793
$measure
NA
odds ratio with 95% C.I. estimate lower upper
Capacete + 1.0000000 NA NA
Capacete - 0.4817645 0.2630993 0.8357342
$p.value
NA
two-sided midp.exact fisher.exact chi.square
Capacete + NA NA NA
Capacete - 0.005088331 0.005688151 0.006857643
$correction
[1] FALSE
attr(,"method")
[1] "Conditional MLE & exact CI from 'fisher.test'"Há mais três métodos disponíveis, de acordo com a documentação da função.
A library epitools também permite o cálculo de \(RR\) com intervalo de confiança. Ainda que não seja adequada para o presente exemplo, verifiquem sua sintaxe. Há três métodos disponíveis, como por exemplo:
library(epitools)
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
epitools::riskratio(tabela, method="wald")$data
Trauma + Trauma - Total
Capacete + 17 138 155
Capacete - 130 508 638
Total 147 646 793
$measure
NA
risk ratio with 95% C.I. estimate lower upper
Capacete + 1.0000000 NA NA
Capacete - 0.8943256 0.8357184 0.9570428
$p.value
NA
two-sided midp.exact fisher.exact chi.square
Capacete + NA NA NA
Capacete - 0.005088331 0.005688151 0.006857643
$correction
[1] FALSE
attr(,"method")
[1] "Unconditional MLE & normal approximation (Wald) CI"e também uma variante com bootstrapping:
library(epitools)
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
epitools::riskratio(tabela, replicates=1e6)$data
Trauma + Trauma - Total
Capacete + 17 138 155
Capacete - 130 508 638
Total 147 646 793
$measure
NA
risk ratio with 95% C.I. estimate lower upper
Capacete + 1.0000000 NA NA
Capacete - 0.8943256 0.8357184 0.9570428
$p.value
NA
two-sided midp.exact fisher.exact chi.square
Capacete + NA NA NA
Capacete - 0.005088331 0.005688151 0.006857643
$correction
[1] FALSE
attr(,"method")
[1] "Unconditional MLE & normal approximation (Wald) CI"| Lembre-se de nunca decidir a respeito de risk ratio ou odds ratio sem considerar o intervalo de confiança. Por isso, evite fazer o cálculo manualmente: use R! |
Dois códigos úteis, com vários dos elementos apresentados neste texto
|
4 Simulação e Gráficos |
1. tabelacontingencia2x2_capacetetrauma.R
Em tabelacontingencia2x2_capacetetrauma.R, que você pode adaptar para as suas necessidades, processamos o mesmo exemplo (procure a função sink() que desvia o resultado da tela para um arquivo texto). Além da estatística \(X^2\), computamos outros valores interessantes. Observe sua saída (tabelacontingencia2x2_capacetetrauma.txt):
source("tabelacontingencia2x2_capacetetrauma.R")
Resultados armazenados em tabelacontingencia2x2_capacetetrauma.txtEm Qui2.R também processamos o mesmo exemplo. Esta versão aplica bootstrapping de um jeito mais manual e produz gráficos associados aos conceitos que discutimos.
Execute com:
source("Qui2.R")
**** TABELA DE CONTINGÊNCIA ****
Considere:
Coluna+ Coluna-
Linha + a b Total lin. +
Linha - c d Total lin. -
Total col. + Total col. - Total geral
Informe:
Nome nas linhas (ate 10 letras): Capacete
Nome nas colunas (ate 10 letras): Trauma
Definido:
Trauma + Trauma -
Capacete + a b Total de Trauma +
Capacete - c d Total de Trauma -
Total de Capacete + Total de Capacete - Total geral
Forneca os valores de a, b, c, d:
a: 17
b: 138
c: 130
d: 508
Defina alfa = prob. erro do tipo I).
(número entre 0 e 1, default = 0.05).
alfa: 0.05
Quantas iteracoes para simular?
(entre um numero inteiro; default n=1e4)
iteracoes: 1e4
Exibir tabelas? (exibir lentifica a simulacao)
0=nao, 1=sim; default eh 0: 0
Observe a simulação em andamento na área de Plots. Quando terminar, procure os gráficos finalizados. Os principais são:
- as distribuições de \(\chi^2\) para a hipótese nula e para o valor observado de \(H_1\):
| O valor de \(\beta\) e o poder associados, aqui, são os valores a posteriori, portanto inúteis para qualquer decisão |
|
Flashback: Recorde a tomada de decisão com a distribuição binomial: A estatística e o formato das distribuições mudam (em função, também, do tipo de variável), mas as áreas de \(\alpha\), \(\beta\) e o raciocínio são os mesmos.
|
- a distribuição, com intervalo de confiança, do odds ratio simulado (compare com o obtido acima), indicando a rejeição de \(H_0\) (o valor 1, de referência, está fora do intervalo) e que seu efeito está entre mínimo a intermediário:
|
Odds ratio é uma medida da intensidade de associação entre duas variáveis nominais dicotômicas de exposição e desfecho. De acordo com SHARPE, D. (2015) Your chi-square test is statistically significant: now what? Practical Assessment, Research & Evaluation, 20(8). É, também, uma medida de tamanho de efeito (não depende do tamanho da amostra, é adimensional e adirecional) e sua intensidade pode ser classificada em:
|
- a distribuição do tamanho de efeito (V de Cramer), mostrando que o efeito, apesar de significante, está entre mínimo e pequeno:
|
O programa já indica o tamanho de efeito de acordo com
Tamanho de efeito:
|
O V de Cramer é uma medida do grau de correlação absoluta entre duas variáveis nominais:- novamente a distribuição de \(\chi^2\), mas aqui exibe somente \(H_0\) e o valor-p, indicando a rejeição de \(H_0\):
2. tabelas LxC
O teste Qui-quadrado opera em tabelas maiores, com \(L\) linhas e \(C\) colunas.
Por exemplo, três quimioterápicos foram comparados quanto ao efeito colateral (náusea). Considere \(\alpha=0.05\). Há diferença entre as drogas?
Os dados obtidos foram:
Nausea Sem_nausea
Droga A 3 5
Droga B 7 2
Droga C 6 3
Em Qui2_LxC.R implementamos este exemplo (versão robusta):
source("Qui2_LxC.R") ComNausea SemNausea
A 3 5
B 7 2
C 6 3
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tcrxc
X-squared = 3.0559, df = NA, p-value = 0.2744
V de Cramer = 0.3428334
Observe que não rejeita-se \(H_0\). Os dois critérios são sempre equivalentes:
- o \(X^2\) calculado é menor que \(\chi^2_c\)
- \(p > \alpha\)
Abra o código para ver como foi implementado este teste.
Comentários sobre os métodos tradicionais
|
5 Condições para o teste Qui-quadrado (métodos não robustos) |
Antes das implementações computacionais e a facilidade do uso do R, a partir de uma tabela 2x2 era necessário construir a tabela dos valores esperados.
- Valores observados:
Trauma + Trauma - Capacete + 17 138 155 Capacete - 130 508 638 147 646 793 - Valores esperados:
Trauma + Trauma - Capacete + 28.7 126.3 Capacete - 118.3 519.7Cada valor esperado é obtido pelo produto das marginais, dividido pelo total geral:
- \((147 \cdot 155) / 793 \approx 28.7\)
- \((147 \cdot 638) / 793 \approx 118.3\)
- \((646 \cdot 155) / 793 \approx 126.3\)
- \((646 \cdot 638) / 793 \approx 519.7\)
Avaliando os dados observados e os valores esperados, era necessário verificar se o teste Qui-quadrado podia ser aplicado. Eram as seguintes condições:
SIEGEL, S. & CASTELLAN Jr., N.J. (1988) Nonparametric statistics for the behavioral sciences. 2ª ed. NY: McGraw-Hill, p. 123:
- O teste exato de Fisher para testar independência é adequado apenas para tabela de contingência 2×2 se N < 20
- Se N entre 20 e 40, o teste Qui-quadrado de Pearson pode ser usado se todas as frequências esperadas são maiores que 5.
- Se N > 40, o teste Qui-quadrado de Pearson com correção de Yates pode ser usado.
COCHAN, W.G. (1954) Some Methods for Strengthening the Common \(\chi^2\) Tests. Biometrics: 10(4): 417-451:
- Um valor esperado mínimo de 1 em alguma célula é permitido, desde que não mais que 20% das células tenham valor abaixo de 5
Então, o valor de \(X^2\) era dado por
\[X^2 = \sum{ ({\text{observado}-\text{esperado}})^2 \over {\text{observado}} }\]
neste exemplo calculado por:
\(X^2 = { {{(17-28.7)^2} \over {28.7}} + {{(138-126.3)^2} \over {126.3}} + {{(130-118.3)^2} \over {118.3}} + {{(508-519.7)^2} \over {519.7}} } \approx 7.31\)
ou, em tabelas 2x2, era recomendado usar a correção de continuidade de Yates, dada por:
\[X^2 = \sum{ (|{\text{observado}-\text{esperado}}|-0.5)^2 \over {\text{observado}} }\]
Neste exemplo:
\(X^2 = { {{(|17-28.7|-0.5)^2} \over {28.7}} + {{(|138-126.3|-0.5)^2} \over {126.3}} + {{(|130-118.3|-0.5)^2} \over {118.3}} + {{(|508-519.7|-0.5)^2} \over {519.7}} } \approx 6.7\)
Então, era necessário consultar uma tabela com os valores críticos de \(\chi^2_c\) previamente calculados para a tomada de decisão (sem obter, portanto, o valor-\(p\)). A tabela permitia escolher o valor de \(\alpha\) e, sabendo-se os graus de liberdade (neste caso \(\nu = 1\)), localizar \(\chi^2_c\) (neste exmeplo igual a \(3.841\)): 
Caso \(\chi^2\) não atendesse as exigências para ser aplicado, a alternativa (se fosse tabela 2x2) era o Teste Exato de Fisher. Este envolve calcular o valor-\(p\) por fatoriais de tabelas progressivamente mais extremas. Por exemplo: 
Neste caso o valor exato de \(p\) era calculado e comparado diretamente com \(\alpha\) para a tomada de decisão.
O que diz a literatura mais recente sobre estes métodos tradicionais e não robustos?
Métodos avançados
Análise post hoc
Esta análise permite, quando se rejeita \(H_0\), localizar quais células da tabela de contingência mais contribuíram para esta rejeição.
exemplo: etilismo e tabagismo
Para descobrir se existe uma relação entre tabagismo e etilismo a partir de um estudo realizado com 100 estudantes universitários, encontrou-se:
| Tabagista | Não tabagista | |
|---|---|---|
| Etilista | 50 | 15 |
| Não etilista | 20 | 25 |
|
A manipulação de dados em R é poderosa, mas nem sempre fácil. Para a entrada de dados na forma de uma matriz 2x2 o seguinte código funciona: Obtendo-se: O problema está na entrada de dados. Somente espaços em branco podem ser usados para a montagem da tabela. O alinhamento é visual, mas não necessário para o R. Funciona igualmente Portanto, pode ser ruim a visualização de sua entrada. Uma opção mais “limpa” é sempre colocar os dados em planilhas no formato .xls ou .xlsx e usar a função readxl::read_excel() para utilizá-los em R: A função read_excel() avisa que teve que suprir um nome (por causa da célula A1 que está vazia), mas esta célula será desprezada adiante. No entanto, temos um problema a resolver. A função read_excel() produz um data frame e a função chisq.test() precisa de uma matriz para sua entrada: Temos, portanto, que converter TCnew para a representação em matrix. Existe uma função para isto: Ainda não é o desejado. TCmatrix tem 2 linhas mas 3 colunas, além de ter perdido os nomes das linhas que TCnew trazia em sua primeira coluna: A solução, portanto, é transferir o conteúdo de TCnew como nomes das linhas de TCmatrix e, então, deletar a coluna 1 de TCmatrix que contém valores NA. Colocando tudo junto, para terminar com uma matriz chamada TC (como era originalmente), podemos ler a planilha em uma variável (e.g., TCtmp), reorganizar seu conteúdo em TC e remover TCtmp da memória com a função remove(). O código enxuto é: Uma mensagem:
apareceu aqui. Warnings são avisos do R que não impedem o funcionamento do programa (avalie, caso a caso, se precisa corrigir algo; aqui não é necessário). Adiante mostramos como é removê-la. |
Esta planilha é Neste exemplo, o objetivo é executar o teste Qui-quadrado mas, no caso de rejeitar-se \(H_0\), executar uma análise post hoc. Observe o código de TabelaContingencia2x2_TabagismoEtilismo.R que pega os dados da planilha tabagismo_e_etilismo_2x2.xlsx:
# TabelaContingencia2x2_TabagismoEtilismo.R
# desabilita warnings
options(warn=-1)
# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- read_excel("tabagismo_e_etilismo_2x2.xlsx")
TC <- data.matrix(TCtmp)
rownames(TC) <-as.character(unlist(TCtmp[1:2,1]))
TC <- TC[,-1]
remove(TCtmp)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=1e6))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("\nResiduos ajustados standardizados corrigidos por momento (MCSTARs):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
print(MCSTAR)
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
cat("\n|MCSTAR critico| (alfaBonferroni=5%/",df,") = ",zcrit,"\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV,
"de dependencia entre as duas variaveis nominais\n")
cat("\nTeste de razao de chances (OR) robusto:\n")
resft <- fisher.test(TC,conf.level = 1-alfa) # Teste de OR robusto
print(resft)
# reabilita warnings
options(warn=0)que produz a seguinte saída:
New names:
* `` -> ...1
------------------------------------------
Dados:
------------------------------------------
Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: TC
X-squared = 12.121, df = NA, p-value = 0.000666
X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 1
Heuristica de significancia (rej. H0 se X^2/gl > 2): 12.12149
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standardizados corrigidos por momento (MCSTARs):
Tabagista NaoTabagista
Etilista 6.963186 -6.963186
NaoEtilista -6.963186 6.963186
|MCSTAR critico| (alfaBonferroni=5%/1) = 1.959964
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.3319569
Grau intermediario de dependencia entre as duas variaveis nominais
Teste de razao de chances (OR) robusto:
Fisher's Exact Test for Count Data
data: TC
p-value = 0.0006368
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.694576 10.331999
sample estimates:
odds ratio
4.107342 A novidade está na sessão Residuos ajustados standardizados corrigidos por momento (MCSTARs). Analise pelos valores positivos acima do valor assinalado como [MCSTAR critico]: são escores \(z\) relativos a quanto cada célula contribuiu para o qui-quadrado estimado. Neste exemplo, a diagonal principal é a responsável pela rejeição de \(H_0\), indicando associação concordante entre etilismo e tabagismo.
|
A única “sujeira” deste código é:
Isto desaparece se colocarmos algo, que será desprezado, na célula A1 da planilha. Experimente rodar o código com a planilha tabagismo_e_etilismo_2x2_v2.xlsx. Note, também, o uso das linhas:
|
Qui-quadrado mínimo
Segundo este autor, a correção de Yates não deve ser usada: 
Note a relação entre esta forma de calcular o qui-quadrado e o odds-ratio: 
Qui-quadrado para variáveis ordinais
O teste qui-quadrado tradicional trata os dois fatores com categorias nominais. Há muitos casos, porém, em que estes fatores são ordinais. Embora não seja incorreto utilizar o qui-quadrado (toda variável ordinal é também nominal), alguma informação está sendo perdida.
Existe uma medida, o Gama de Goodman-Kruskal, que considera ambas as variáveis como categóricas ordinais, potencialmente melhorando o resultado do teste. É uma estatistica nao-parametrica que mede a força de associação em dados cruzados em tabela, com as hipóteses:
\[H_o: \gamma = 0\] \[H_o: \gamma \ne 0\]
Neste estudo, Goodman (1979) desenvolve a análise com variáveis ordinais. 
Os dados desta tabela são, originalmente, de:
Duncan, O.D. (1979), How Destination Depends on Origin in the Occupational Mobility Table. American Journal of Sociology, 84, 793-803.
As categorias são as seguintes:
- Professional and high administrative
- Managerial and executive
- Inspectional, supervisory, and other nonmanuar (high grade)
- Inspectional, supervisory, and other nonmanual (lower grade)
- (V)(a) Routine grades of nonmanual
- (V)(b) Skilled manual
- Semiskilled manual
- Unskilled manual
Portanto, a variável de status ocupacional é, claramente, qualitativa ordinal.
| O repositório do R tem mais que pacotes. Há vários bancos de dados armazenados nele, para uso em exemplos. Aqui usaremos os dados deste estudo, disponíveis no pacote “dataset”. |
Implementado em TCrXc_StatusOcupacional.R:
# TCrXc_StatusOcupacional.R
library(datasets)
library(DescTools)
eiras.show.MCSTAR <- function (MCSTAR, zcrit)
{
MCSTARtxt <- round(MCSTAR,3)
for (r in 1:nrow(MCSTAR))
{
for (c in 1:nrow(MCSTAR))
{
if (MCSTAR[r,c] > zcrit)
{
MCSTARtxt[r,c] <- paste("#",MCSTARtxt[r,c],sep="")
} else
{
MCSTARtxt[r,c] <- paste(" ",MCSTARtxt[r,c],sep="")
}
}
}
return(MCSTARtxt)
}
# desabilita warnings
options(warn=-1)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
# Duncan, O.D. (1979), How Destination Depends on Origin in the
# Occupational Mobility Table. American Journal of Sociology 84, 793-803.
tcrxc <- datasets::occupationalStatus
# para recolocar as categorias originais
categories <- c("I","II","III","IV","Va","Vb","VI","VII")
rownames(tcrxc) <- categories
colnames(tcrxc) <- categories
prmatrix(tcrxc)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(tcrxc,simulate.p.value=TRUE,B=1e6))
N <- sum(res$observed)
nL <- nrow(tcrxc) # ou dim(TC)[1]
nC <- ncol(tcrxc) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=df), "\n", sep="")
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("\nResiduos ajustados standarizados corrigidos por momento (MCSTARs):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
print(eiras.show.MCSTAR(MCSTAR,zcrit))
cat("\n|MCSTAR critico| (alfaBonferroni=5%/",df,") = ",zcrit,"\n", sep="")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV,
"de dependencia entre as duas variaveis nominais\n")
# Goodman, L. A. (1979) Simple Models for the Analysis of Association
# in Cross-Classifications having Ordered Categories.
# J. Am. Stat. Assoc., 74 (367), 537–552.
gama <- DescTools::GoodmanKruskalGamma(tcrxc)
cat("\nGama de Goodman-Kruskal = ", gama,"\n", sep="")
# EDWARDES, MD & BALTZAN, M (2000) The generalization of the odds ratio,
# rik ratio and risk difference to rxk tables.
# Statistics in Medicine 19:1901-14.
ORg <- (1+gama)/(1-gama)
cat("\nOR generalizado = (1+",gama,")/(1-",gama,") = ", ORg,"\n", sep="")
# reabilita warnings
options(warn=0)A saída é:
------------------------------------------
Dados:
------------------------------------------
I II III IV Va Vb VI VII
I 50 19 26 8 7 11 6 2
II 16 40 34 18 11 20 8 3
III 12 35 65 66 35 88 23 21
IV 11 20 58 110 40 183 64 32
Va 2 8 12 23 25 46 28 12
Vb 12 28 102 162 90 554 230 177
VI 0 6 19 40 21 158 143 71
VII 0 3 14 32 15 126 91 106
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tcrxc
X-squared = 1416, df = NA, p-value = 1e-06
X^2 critico de 95% = 66.33865
Graus de liberdade (nao fornecidos por bootstrapping): 49
Heuristica de significancia (rej. H0 se X^2/gl > 2): 28.89877
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standarizados corrigidos por momento (MCSTARs):
destination
origin I II III IV Va Vb VI VII
I #28.022 #6.466 #4.851 -2.711 -0.804 -7.091 -4.336 -4.284
II #6.535 #15.194 #6.477 -0.475 0.201 -6.217 -4.43 -4.437
III 0.706 #6.01 #7.195 #3.979 2.782 -3.966 -6.129 -4.134
IV -1.369 -0.926 1.7 #6.772 0.826 0.847 -3.453 -5.132
Va -1.436 0.409 -0.87 0.701 #5.188 -1.363 0.388 -1.982
Vb -6.546 -6.397 -3.505 -1.856 -0.703 #7.926 0.031 1.551
VI -4.57 -4.075 -4.744 -3.41 -2.462 0.329 #9.978 2.718
VII -4.152 -4.315 -4.744 -3.427 -2.901 -0.678 #4.169 #11.153
|MCSTAR critico| (alfaBonferroni=5%/49) = 3.284839
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.2404799
Grau pequeno de dependencia entre as duas variaveis nominais
Gama de Goodman-Kruskal = 0.4209053
OR generalizado = (1+0.4209053)/(1-0.4209053) = 2.453667É praticamente o mesmo código utilizado para o exemplo de tabagismo e etilismo, com algumas diferenças:
- adicionamos uma função, eiras.show.MCSTAR() para exibir a matrix resultante de MCSTAR com menor número de casas decimais e símbolos adicionais (“#”) para que os valores acima de \(z\) crítico sejam mais facilmente localizados.
- em uma tabela com 8 linhas e 8 colunas não existe definição para o odds ratio. No lugar deste, aparece o gama de Godman-Kruskal e um cálculo de OR generalizado.
Apenas para comparação, o teste qui-quadrado comum produziria o seguinte:
source("TCrXc_StatusOcupacional.R")
------------------------------------------
Dados:
------------------------------------------
I II III IV Va Vb VI VII
I 50 19 26 8 7 11 6 2
II 16 40 34 18 11 20 8 3
III 12 35 65 66 35 88 23 21
IV 11 20 58 110 40 183 64 32
Va 2 8 12 23 25 46 28 12
Vb 12 28 102 162 90 554 230 177
VI 0 6 19 40 21 158 143 71
VII 0 3 14 32 15 126 91 106
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tcrxc
X-squared = 1416, df = NA, p-value = 1e-06
X^2 critico de 95% = 66.33865
Graus de liberdade (nao fornecidos por bootstrapping): 49
Heuristica de significancia (rej. H0 se X^2/gl > 2): 28.89877
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standarizados corrigidos por momento (MCSTARs):
destination
origin I II III IV Va Vb VI VII
I #28.022 #6.466 #4.851 -2.711 -0.804 -7.091 -4.336 -4.284
II #6.535 #15.194 #6.477 -0.475 0.201 -6.217 -4.43 -4.437
III 0.706 #6.01 #7.195 #3.979 2.782 -3.966 -6.129 -4.134
IV -1.369 -0.926 1.7 #6.772 0.826 0.847 -3.453 -5.132
Va -1.436 0.409 -0.87 0.701 #5.188 -1.363 0.388 -1.982
Vb -6.546 -6.397 -3.505 -1.856 -0.703 #7.926 0.031 1.551
VI -4.57 -4.075 -4.744 -3.41 -2.462 0.329 #9.978 2.718
VII -4.152 -4.315 -4.744 -3.427 -2.901 -0.678 #4.169 #11.153
|MCSTAR critico| (alfaBonferroni=5%/49) = 3.284839
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.2404799
Grau pequeno de dependencia entre as duas variaveis nominais
Gama de Goodman-Kruskal = 0.4209053
OR generalizado = (1+0.4209053)/(1-0.4209053) = 2.453667| Especificamente para tabelas 2x2 existe o \(Q\) de Yule, dado por \[Q = {{ad-bc} \over {ad+bc}}\] Assim como o Gama de Goodman-Kruskal, varia de \(-1\) a \(+1\), mas como se aplica a tabelas 2x2, relaciona-se com o odds-ratio: \[Q = {{OR-1} \over {OR+1}}\] Uma função está implementada em psych::Yule(x, Y=FALSE), que devolve o \(Q\) de Yule quando \(x\) é uma tabela de frequências absolutas 2x2. |
Qui-quadrado para tendência de Cochran-Armitage
Para situações em que uma variável é ordinal e a outra é nominal, por exemplo quando há um modelo de dose-resposta com um desfecho dicotômico associado a uma exposição que tem mais que dois níveis, é preferível utilizar o teste conhecido como Cochran-Armitage test for trend.
|
|
exemplo: espessura da dobra cutânea e menarca
Os pesquisadores Beckles et al. (1985) estudaram a associação entre a variável de a faixa etária da menarca (\(<\) 12 anos e $$12 anos) e o fator a espessura da dobra cutânea do tríceps (pequena, intermediária e grande).
|
|
Caso ambas as variáveis sejam consideradas nominais, o teste adequado é o qui-quadrado de Pearson. Caso a variável considerada como exposição seja ordinal e o desfecho seja dicotômico (tabela \(K\)x\(2\)), o teste qui-quadrado de Cochran Armitage pode ser usado.
Para comparação, ambos os testes são executados em Teste qui-quadrado de Cochran-Armitage.R:
# Teste qui-quadrado de Cochran-Armitage.R
# Cochran-Armitage test for trend
# Menarca precoce & espessura da dobra cutânea do tríceps
# Fonte: Beckles et al. (1985) International Journal of Obesity 9:127-35,
# apud Kirkwood & Sterne (2003): Chapter 17: Chi-squared tests for 2x2 and
# larger contingency tables.
library(readxl)
library(DescTools)
# desabilita warnings
options(warn=-1)
# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- read_excel("menarca.xlsx")
TC <- data.matrix(TCtmp)
rownames(TC) <-as.character(unlist(TCtmp[,1]))
TC <- TC[,-1]
remove(TCtmp)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=1e6))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\nTeste qui-quadrado de Cochran-Armitage:\n")
res <- DescTools::CochranArmitageTest(TC)
print(res)
x2ca <- res$statistic^2
p <- res$p.value
cat("\nX-squared trend =",x2ca, ", df = 1, ", "p-value =", p,"\n")
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV,
"de dependencia entre as duas variaveis nominais\n")
# reabilita warnings
options(warn=0)que produz:
------------------------------------------
Dados:
------------------------------------------
menor12 maiorigual12
Pequena 15 156
Intermediaria 29 197
Grande 36 186
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: TC
X-squared = 4.7594, df = NA, p-value = 0.09404
X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 2
Heuristica de significancia (rej. H0 se X^2/gl > 2): 2.379692
Teste qui-quadrado de Cochran-Armitage:
Cochran-Armitage test for trend
data: TC
Z = 2.1783, dim = 3, p-value = 0.02938
alternative hypothesis: two.sided
X-squared trend = 4.744926 , df = 1, p-value = 0.02938481
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.08768596
Grau minimo de dependencia entre as duas variaveis nominaisCompare os valores \(p\) de ambos os testes. Cochran-Armitage, ao considerar a espessura da prega cutânea como “dose”, ordinal, indica significância estatística onde o teste convencional não detectava.
Cochran-Armitage generalizado
Para desfecho com três ou mais categorias (os desfechos são tratados como variáveis nominais) existe função no pacote coin::chisq.test
Um exemplo está em Teste qui-quadrado de Cochran-Armitage generalizado.R:
library(coin)
library(rcompanion)
# The chisq_test function in the coin package can be used to conduct a test of
# association for a contingency table with one ordered nominal variable and
# one non-ordered nominal variable.
# The Cochran–Armitage test is a special case of this where the
# non-ordered variable has only two categories.
# The scores option is used to indicate which variable should be treated as
# ordered, and the spacing of the levels of this variable.
Job <- matrix(c(1,2,1,0, 3,3,6,1, 10,10,14,9, 6,7,12,11), 4, 4,
dimnames = list(income = c("< 15k", "15-25k", "25-40k", "> 40k"),
satisfaction = c("VeryD", "LittleD",
"ModerateS", "VeryS")))
Job <- as.table(Job)
print(coin::chisq_test(Job, scores = list("income"= c(1, 2, 3, 4))))
# Teste post hoc
print(ph <- rcompanion::pairwiseOrdinalIndependence(Job,compare="column"))
try(rcompanion::cldList(p.value ~ Comparison, data = ph, threshold = 0.05))que produz:
Loading required package: survival
Attaching package: 'survival'The following object is masked from 'package:epitools':
ratetable
Asymptotic Generalized Pearson Chi-Squared Test
data: satisfaction by
income (< 15k < 15-25k < 25-40k < > 40k)
chi-squared = 3.1362, df = 3, p-value = 0.3711
Comparison p.value p.adjust
1 VeryD : LittleD 0.451 0.541
2 VeryD : ModerateS 0.351 0.526
3 VeryD : VeryS 0.161 0.526
4 LittleD : ModerateS 0.698 0.698
5 LittleD : VeryS 0.242 0.526
6 ModerateS : VeryS 0.271 0.526Qui-quadrado de Mantel-Haenszel
Este teste utiliza uma tabela tridimensional. Neste exemplo temos a contagem de mortes entre fumantes e não fumantes estratificadas por faixa etária. Esta versão de Qui-quadrado faz o teste global e, depois, estratifica pela terceira variável (faixa etária, neste exemplo).
| Faixa Etária | 18-24 | 25-34 | 35-44 | 45-54 | 55-64 | 65-74 | 75 | |
|---|---|---|---|---|---|---|---|---|
| Tabagista | Desfecho | |||||||
| Sim | Morta | 2 | 3 | 14 | 27 | 51 | 29 | 13 |
| Viva | 53 | 121 | 95 | 103 | 64 | 7 | 0 | |
| Não | Morta | 1 | 5 | 7 | 12 | 40 | 101 | 64 |
| Viva | 61 | 152 | 114 | 66 | 81 | 28 | 0 |
Os dados precisam ser reorganizados para que construamos uma tabela tridimensional requerida pela funções relacionadas a este teste. Veja como os dados foram acomodados na planila fumo_e_faixaetaria.xlsx.
Implementamos o teste em Teste Qui-quadrado de Mantel-Haenszel.R:
# Teste Qui-quadrado de Mantel-Haenszel.R
# desabilita warnings
options(warn=-1)
# Tabela de contingencia 2x2 segmentada
# Teste robusto da razao de chances (OR) de Mantel-Haenszel
# Relacao entre tabagismo e sobrevivência em 20 anos (1974-1994)
# em 1.134 mulheres adultas do Reino Unido
# Delineamento: coorte
# Fonte: APPLETON, D. R. et al. (1996) Ignoring a covariate:
# An example of Simpson's paradox. The American Statistician,
# 50(4): 340-1.
# A aspa inicial TEM que comecar na primeira coluna da linha
# e os espacamentos distintos dos dessa tabela podem causar
# problemas de na geracao da tabela horizontalizada (ftable).
library(vcd)
library(rcompanion)
library(coin)
library(readxl)
library(tidyverse)
TCtmp <- read_excel("fumo_e_faixaetaria.xlsx")
TCtmp <- gather(TCtmp, var, Contagem, -Idade)
TCtmp <- separate(TCtmp, var, c('Exposicao','Desfecho'))
TCS <- TCtmp %>% xtabs(Contagem ~ Exposicao + Desfecho + Idade, .)
print(TCS)
# Mantel-Haenszel test of the null that two nominal variables
# are conditionally independent in each stratum,
# assuming that there is no three-way interaction.
# Woolf test for homogeneity of ORs across strata on 2 x 2 x k tables:
# If significant, M-H test is not appropriate.
print(ftable(TCS)) # Display a flattened table (tabela horizontalizada)
cat("\nHomogeneidade dos ORs:\n")
print(vcd::woolf_test(TCS))
# Teste para tabela 2x2xK
cat("\nQui-quadrado de Mantel-Haenszel:\n")
print(mantelhaen.test(TCS, exact=TRUE)) # Teste exato de OR de Mantel-Haenszel
# Teste para tabela LxCxK
cat("\nTeste de Cochran-Mantel-Haenszel (coin::cmh_test):\n")
print(coin::cmh_test(TCS, distribution = approximate(nresample = 1e6)))
cat("\nAnalise estratificada por idade:\n")
print(rcompanion::groupwiseCMH(TCS,
group = 3,
fisher = TRUE,
gtest = FALSE,
chisq = FALSE,
method = "fdr",
correct = "none",
digits = 3))
cat("\n")
print(ors <- vcd::oddsratio(TCS, log=FALSE)) # Show OR for each 2x2
cat("\n")
lnors <- vcd::oddsratio(TCS, log=TRUE) # Show ln(OR) for each 2x2
print(lnors)
print(summary(lnors))
print(confint(lnors))
plot(lnors, xlab = "Faixa Etaria", ylab = "ln(OR)")
# reabilita warnings
options(warn=0)que gera a seguinte saída:
Loading required package: grid
Attaching package: 'vcd'The following object is masked from 'package:epitools':
oddsratio── Attaching packages ──────────────────────────────────────────────────── tidyverse 1.3.0 ──✓ ggplot2 3.3.0 ✓ purrr 0.3.3
✓ tibble 3.0.1 ✓ dplyr 0.8.5
✓ tidyr 1.0.2 ✓ stringr 1.4.0
✓ readr 1.3.1 ✓ forcats 0.5.0── Conflicts ─────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag(), , Idade = 18-24
Desfecho
Exposicao Morta Viva
NaoTabagista 1 61
Tabagista 2 53
, , Idade = 25-34
Desfecho
Exposicao Morta Viva
NaoTabagista 5 152
Tabagista 3 121
, , Idade = 35-44
Desfecho
Exposicao Morta Viva
NaoTabagista 7 114
Tabagista 14 95
, , Idade = 45-54
Desfecho
Exposicao Morta Viva
NaoTabagista 12 66
Tabagista 27 103
, , Idade = 55-64
Desfecho
Exposicao Morta Viva
NaoTabagista 40 81
Tabagista 51 64
, , Idade = 65-74
Desfecho
Exposicao Morta Viva
NaoTabagista 101 28
Tabagista 29 7
, , Idade = 75+
Desfecho
Exposicao Morta Viva
NaoTabagista 64 0
Tabagista 13 0
Idade 18-24 25-34 35-44 45-54 55-64 65-74 75+
Exposicao Desfecho
NaoTabagista Morta 1 5 7 12 40 101 64
Viva 61 152 114 66 81 28 0
Tabagista Morta 2 3 14 27 51 29 13
Viva 53 121 95 103 64 7 0
Homogeneidade dos ORs:
Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)
data: TCS
X-squared = 3.2061, df = 6, p-value = 0.7826
Qui-quadrado de Mantel-Haenszel:
Exact conditional test of independence in 2 x 2 x k tables
data: TCS
S = 230, p-value = 0.01591
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
0.4538410 0.9355508
sample estimates:
common odds ratio
0.6534856
Teste de Cochran-Mantel-Haenszel (coin::cmh_test):
Approximative Generalized Cochran-Mantel-Haenszel Test
data: Desfecho by
Exposicao (NaoTabagista, Tabagista)
stratified by Idade
chi-squared = 5.8443, p-value = 0.01584
Analise estratificada por idade:
Group Test p.value adj.p
1 18-24 Fisher 0.6000 1.000
2 25-34 Fisher 1.0000 1.000
3 35-44 Fisher 0.0706 0.291
4 45-54 Fisher 0.3650 0.852
5 55-64 Fisher 0.0830 0.291
6 65-74 Fisher 1.0000 1.000
7 75+ Fisher 1.0000 1.000
odds ratios for Exposicao and Desfecho by Idade
18-24 25-34 35-44 45-54 55-64 65-74 75+
0.5219512 1.2519906 0.4314109 0.7074504 0.6223718 0.9054416 4.7777778
log odds ratios for Exposicao and Desfecho by Idade
18-24 25-34 35-44 45-54 55-64 65-74
-0.65018114 0.22473479 -0.84069420 -0.34608770 -0.47421763 -0.09933253
75+
1.56397554
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
18-24 -0.650181 1.049580 -0.6195 0.53561
25-34 0.224735 0.694493 0.3236 0.74624
35-44 -0.840694 0.470642 -1.7863 0.07406 .
45-54 -0.346088 0.375584 -0.9215 0.35681
55-64 -0.474218 0.268109 -1.7687 0.07694 .
65-74 -0.099333 0.460621 -0.2156 0.82926
75+ 1.563976 2.022270 0.7734 0.43930
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2.5 % 97.5 %
18-24 -2.7073204 1.40695808
25-34 -1.1364462 1.58591573
35-44 -1.7631351 0.08174672
45-54 -1.0822181 0.39004270
55-64 -0.9997024 0.05126713
65-74 -1.0021328 0.80346776
75+ -2.3996018 5.52755287